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Abstract. The objective of the paper is to evaluate the ability of reactive transport models and their nu- 
merical implementations (such as MIN3P) to simulate simple microbial transformations in conditions of 
chemostat or gradostat models, that are popular in microbial ecology and waste treatment ecosystems. To 
make this comparison, we first consider an abstract ecosystem composed of a single limiting resource and 
a single microbial species that are carried by advection. In a second stage, we consider another microbial 
species in competition for the same limiting resource. Comparing the numerical solutions of the two models, 
we found that the numerical accuracy of simulations of advective transport models performed with MIN3P 
depends on the evolution of the concentrations of the microbial species: when the state of the system is close 
to a non-hyperbolic equilibrium, we observe a numerical inaccuracy that may be due to the discretization 
method used in numerical approximations of reactive transport equations. Therefore, one has to be cautious 
about the predictions given by the models. 

Key-words. : reactive transport models, chemostat model, microbial growth, numerical simulation. 



1 Introduction 

The chemostat is a popular apparatus, invented simultaneously by Monod [29] and Novick and Szilard 32 , 
that allows the continuous culture of micro-organisms in a controlled medium. The chemostat has the ad- 
vantage to study bacteria growth at steady state, in contrast to batch cultivation. The chemostat model 
serves also as a representation of aquatic natural ecosystems such as lakes. In the classical experiments 
involving chemostats, the medium is assumed to be perfectly mixed, that justifies mathematical models de- 
scribed by systems of ordinary differential equations [36| . In natural ecosystems, ground- waters or industrial 
applications that use large reservoirs, the assumption of perfectly mixed medium is questionable, leading 
to spatialized models such as systems of nonlinear partial differential equations |37| . However, nonlinear 
partial differential equations are difficult mathematical objects to understand, analyze and simulate. Even 
numerical schemes pose significant difficulties, particularly when solving coupled systems involving stiff re- 
actions [3l |33l 04] . Spatial considerations can be introduced in the "classical" model of the chemostat in 
simpler ways, as it is done in the gradostat model [18j that mimics a series of interconnected chemostats of 
identical volumes. Other kinds of interconnection can be considered in order to cope with heterogeneity of 
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porous media, considering stagnant water with small diffusion, mixing due to porosity [131 [301 - ^ course, 
such representations are still oversimplified with regard to the complexity of natural porous media. 

Over the past three decades, numerous reactive transport codes have been developed to study complex 
interactions between geochemical and transport processes in porous media. A number of reactive transport 
computer codes exist. Let us mention for instance the models COMEDIE-2D [6], CRUNCH [35], PHREEQC 
[3J, ECOSAT [15], ORCHESTRA [25] , RAFT g], RT3D [5J, HYTEC 02], HP1 Q3] and MIN3P [22] . To 
our knowledge, some of these numerical tools, such as COMEDIE-2D and PHREEQC, are not suitable 
for unsaturated porous media and thus cannot be readily applied to soils (excepted in the special case of 
wetland soils). A range of other limitations can be found as well. For example, RT3D does not include 
equilibrium-controlled reactions, while ECOSAT neglects kinctically-controlled reactions and is limited to a 
single spatial dimension. 

The standard approach for evaluating the computational accuracy of a reactive transport code is to 
compare its numerical results to those obtained from an analytical or a semi-analytical solution [9l l33l l4ll [43] . 
Unfortunately, analytical solutions are only available for simplified systems, such as the reactive transport of 
a single solute in 1-D homogeneous systems at steady state, which is well behind the actual capabilities of the 
models. To remediate to this, inter-comparison of numerical codes has been largely employed. This inter- 
comparison involves the independent solution of the same problem using a variety of numerical techniques 

En n ma es El- 

This study aims at comparing the accuracy of a reactive transport model with other kinds of models such 
as the mathematical model of the chemostat. This confrontation takes place in the framework of microbial 
ecology, for which concepts of competition and coexistence are crucial [TJ [2J [5J [35J [321 E0] . We have chosen the 
reactive transport code MIN3P [22, for this study. This reactive transport model is notably recognized for its 
numerical robustness [H[22]. In addition, the model MIN3P can simulate general reactive transport problems 
in variably saturated media for ID to 3D systems. The flow solution is based on Richard's equation [31], and 
solute transport is simulated by means of the advection-dispersion equation. Gas transport can be taken into 
account as well, either by considering advection and Fick diffusion or the Dusty Gas Model [251127]. A range of 
bio-geochemical processes are included in MIN3P (aqueous speciation, mineral dissolution-precipitation, gas 
dissolution/exsolution, ion exchange, and competitive or non competitive sorption). A generalized kinetic 
expression for dissolution-precipitation and intra-aqueous reactions allows the consideration of fractional 
order or Monod-type rate expressions, and parallel reaction pathways. This code has been used for a number 
of applications in different fields of environmental science, ranging from inorganic and organic contaminant 
transport and groundwater remediation [201 O] H31 [25] to soil hydrology and bio-geochemical cycles in 
terrestrial ecosystems [TTJ [T2J [TpJ . 

The chemostat model is the simplest mathematical model for describing the dynamics of microbial growth 
under a constant flow of substrate, and its theory is well understood [37]. In this work, we consider a series 
of interconnected chemostats for the simulation a one-dimensional heterogeneity, that we compare with the 
solutions provided by reactive transport models considering the same spatial structure. 

2 Material and method 

In practice, a chemostat in laboratory is an apparatus that consists of three connected vessels as shown in 
Fig. [T] The leftmost vessel is called the feed bottle and contains all of the nutrients needed for the growth 
of a microorganism. The central vessel is called the culture vessel, while the last is the overflow or collection 
vessel. The content of the feed bottle is pumped at a constant rate into the culture vessel, while the content 
of the culture vessel is pumped at the same constant rate into the collection vessel. We denote by Si n [mol/l] 
the constant concentration of nutrient pumped with a volumetric flow rate that we denote by Q\t.s" ]. The 
dilution rate [s _1 ] is defined as D = Q/V, where V is the volume of the culture vessel. We shall also denote 
by /u(.)[s _1 ] the specific growth rate of micro-organisms and k the yield factor of the bio-conversion. The 
dynamical model of the chemostat can then been written in the following way. 
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Because of the boundary conditions (i.e. input of nutrients in the culture vessel and output of contents 
from the culture vessel) , numerical implementations of reactive transport models such as MIN3P are not able 
to simulate straightforwardly a single ecosystem such as the one in the culture vessel. Indeed, the use of the 
logarithm of the concentrations in the numerical code prevents to having a null concentration of biomass at 
the input of the culture vessel. In order to simulate an ecosystem in a single tank, one has to consider three 
control volumes. Moreover, for intrinsic reasons, three control volumes is the minimum number for a one 
dimensional discretization in numerical implementations of reactive transport models, such as the MIN3P 
code that we have chosen to simulate our ecosystem. In this way, the numerical implementation is closer 
from the true laboratory experiment with three vessels, that we described above. Nevertheless, we shall refer 
in the following to the chemostat for the culture vessel only. 



Figure 1: A schematic view of the chemostat experiment 



-Bi + Di(Si-i — Si) 



With MIN3P we began by simulating a simple example of three chemostats in series having the same 
volume, in presence of a single biomass B and a single substrate S with the same initial conditions in the 
three control volumes. The dynamical model representing n chemostats connected in series is given by the 
equations: 

dSj _ _a*(Sj) 
~dt ~ 

d ^ = { l i{S i )-D i )B i + D i B i _ u 

where Di = Q/Vi, Si(respBi) represents the substrate concentration (resp biomass concentration) in the ith 
bioreactor (i = 1, , n). So = Si n and Bo = and Vj = — . 

We consider that the qualitative behavior of this system of ordinary differential equations (ODE) is today 
well understood (see the Appendix), and we used different robust numerical schemes for solving this system 
of ODE (Rungc-Kutta, LSODA ...) that were all thoroughly consistent with the analytic analysis of the 
steady states and their stability. That amounts to assume that we can trust the numerical solutions obtained 
by the numerical integration of ODEs for this model. 



For the simulation of the chemostat model with reactive transport models, we consider a boundary 
domain in three dimensions. The boundary conditions for the liquid flow are of first type with a value of 
in the output face, and of second type specifying a flux of Q[l.s~ ] in the input and output faces. A specific 
choice of flow condition gave us a fully water saturated medium at any time . The boundary conditions for 
the reactive part are of second type in the output face and of third type in the input face with a value of 
the substrate concentration equal to 5 in . 

To simulate several chemostats connected in series without diffusion, we have chosen the free diffusion 
coefficient in water and air, as well as the specific storage coefficient equal to zero. The porosity of the 
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Figure 2: Configuration of three control volumes 



medium formed by one domain only is chosen equal to one. The day has been chosen as the time unit, with 
a maximum time increment of 10 -3 day and a starting time step of 10 -10 day. 

Finally, we have introduced a simple theoretical irreversible reaction expressing the bio-conversion of 
moles of substrate into one mole of biomass. The specific growth rate of biomass considered here is of 
Monod's type: 

where /i ma x[s _1 ] represents the maximum of the intra-aqueous kinetic reaction and k s [mol/l] the half satu- 
ration constant. 

In a second stage, we consider the same spatial considerations but with two species instead of one, 
assuming that each species has a growth function that follows a Monod law, as described above. We assume 
that their interaction is due only to a common limiting resource, that is the species compete for the same 
substrate. We focus on the case of a "true" competition, in the sense that we assume that a species is the 
most efficient one when the resource is very rare, while the other species is better when the resource is less 
rare. 



3 Results and discussion 

Denote by S* M (resp. Sq) the value of the substrate concentration computed by MIN3P (resp. by the chemo- 
stat model) at the steady state, and by B* M the value of the biomass concentration at steady state given 
by MIN3P. Let 5 be the absolute value of the difference between and Sq, that serves as an indicator of 
divergence between simulations of both models. 

For the comparison of the two models, we first represent the indicator S as a function of Q. For the study 
of the effect of a spatial discretization, we shall then consider 5 as a function of the number of cells. 
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For convenience, we shall omit the units of the numerical values that we give below (concentrations are 
assumed to be measured in mol.l , volumes in I, flux in Z.s -1 and growth rates in s _1 . 

3.1 The single species case 

Choosing fi max = 4.10~ 5 , we have studied the variation of S with respect to Q in the first and third control 
volumes. Each volume has been chosen equal to 1. On the graphs of 6 that correspond to the first (resp. 
third) control volume (see Figure [3j Output (resp. Input)), we notice that for 10~ 6 < Q < 10~ 5 , one faces 
a significant difference between simulations of both models in the input control volume, and that is even 
greater in the output control volume. But for 1CP 5 < Q < 15. 1CP 6 , one has almost no difference in the 
input control volume, and one can observe a jump of 6 about Q = 10~ 5 in the output control volume. For 
the simulations, we have chosen Si n — 3, k s = 1, B(0) = 2 and 5(0) = 5. 



We know from the theory of the chemostat (see Appendix, Proposition 4.4 1, that for a "perfectly mixed" 
tank, with a single species growing on a single limiting substrate, the condition A '^"- ) > I ensures that the 
biomass B is not wash-out. This result is surprisingly not observed in simulations with MIN3P. To show 
that, we studied the variation of the biomass concentration B with respect to Q in the input control volume. 
We notice that the biomass is washed out for Q > 8.10~ 6 (see Figure [3j B* M in dashed line). But when 
3.10~ 6 < Q < 10~ 5 , we have M< ^'") = 42-— > 1, and one can observe on the graph of S and the difference 



between MIN3P and the chemostat simulations in the input control volume. 

For the case of three chemostats connected in series with the same volume and traversed by the same 
flow rate Q, the removal of the biomass in the input control volume has to lead theoretically to its removal in 
the output control volume. But this is not the case with MIN3P and we can notice on Fig. [3] (Output), that 
for 8.10 -6 < Q < 15. 10 -6 , we obtain the wash-out of biomass in the input control volume, the biomass in 
the output control volume being not yet washed out. In other words, under certain conditions, the microbial 
growths predicted by both simulations are radically different. 



B 



Output 




Figure 3: Comparison for three control volumes 



For the study of the effect of a spatial discretization on the difference between numerical reactive trans- 
port and chemostat models at steady state, we took the same conditions as before, with a maximal kinetic 
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rate equal to fi max = 5.1CP 4 and a flow Q = 1CP 5 . We did vary the number of discretization steps between 
n — 3 and n — 50 and studied the variation of S with respect to n (the total number of cells) in the input 
and output control volumes. Denote by D n the dilution rate in each control volume for a discretization in 
n cells. On the graphs of <5 that corresponds to the input (resp. output) control volume (see Fig. |4j Input 
(resp. Output)), we notice for 3 < n < 37 (that theoretically corresponds to the survival of biomass, because 
one has = H > 1) that we have the same behavior of 5 as before. Similarly for 39 < n < 50 (values 

that correspond to < 1), we have no difference between MIN3P and the chemostat model in the input 

and output control volumes. But for n = 38, we observe the same jump of d in the output control volume 
as previously observed. 

We have observed on this example a significant difference between reactive transport and chemostat 
models when passing from a steady state of survival of the biomass to the wash out steady state, this 
corresponds to a bifurcation passing from two equilibriums (wash-out and biomass survival) to a single 
equilibrium (wash-out). The limiting case corresponds to a non- hyperbolic equilibrium (see the Appendix 
for the definition of hyperbolic equilibriums). Because of the use of the logarithm in the MIN3P code, 
we expect the internal solution to take very large values when the concentration of biomass tends to zero. 
One can also detect this phenomenon on the simulations when noticing a "time lag" between MIN3P and 
chemostat trajectories about the steady state. 

f Output 

30 - ^ D M 




5 10 15 20 25 30 35 40 45 50 

Number of control volume 
Figure 4: Comparison with several control volumes 



3.2 The two species case 

To emphasize the problem that occurs about non-hyperbolic equilibrium in the case of one species, we 
present in this part a more subtle situation, considering two species that compete for the same substrate. 



G 



The extension of the model ([I]) is given by the following equations: 

k\ k 2 

( f i 1 (S)-D)B 1 (2) 
(» 2 (S)-D)B 2 

where functions and n 2 (.) denote the kinetics of species B\ and B 2 respectively. For this system, the 

Proposition |4.5| given in the Appendix shows that non- hyperbolic equilibrium points could exists away from 
the wash-out equilibrium. For the one species case, a non-hyperbolic equilibrium could exists but on the 
boundary of the positive domain only. For the two species case, many non-hyperbolic equilibriums could 
exists on the interior of the positive domain. 



dS 
~dt 
dB l 
~dT 
dB 2 
dt 



Let Xi be the positive solution, when it exists, of ^i(S) = D for i = 1,2. Under the condition Si n > 
max(Ai, A2), the system ^ admits generically three equilibrium points, given by E o (Si n ,0, 0), Ei(Xi, ki(Si n — 
Ai),0) and E 2 (X 2 ,0,k 2 (Si n — X 2 )). We distinguish now two different cases. If for some i = 1,2 one has 



Xi = Sin, then the equilibrium Eo is non- hyperbolic as before (see Proposition 4.4 in the Appendix). More- 
over if one has Ai = A2 < Si n , then E\ = E 2 is also a non hyperbolic equilibrium (Proposition 4.5 in the 
Appendix). Starting from our observations in the case of one species, we have built suitable examples to 
study the behavior of MIN3P about those non-hyperbolic equilibriums. 

To understand the comparison, we first recall that the mathematical theory of the model ^ predicts 
the competitive exclusion in the generic case, that is at most one competitor avoids the extinction (see the 
Appendix). This property refers to the well-known Competitive Exclusion Principle in ecology, that has 
been widely studied in the literature (see for instance [IJ El [39] ) . For the chemostat model, the Principle 
can be stated as follows. 

Considering two increasing growth rates /ii(.) and fJ, 2 (.) such that both Ai and A2 are smaller than Si n 
(for a sufficiently large Si n ). Then, one has the following issue of the competition for large times 

- when Ai < A2, the species B\ avoids the extinction, 

- when Ai > A2, the species B 2 avoids the extinction. 

For the non generic case Ai = A2, it is possible to predict the coexistence of the two species, invalidating the 
Principle (on this single chemostat case). 



In the simulations, we have considered two species B\ and B 2 with a specific growth rate given by 

= 1-10- 3 ^ and M2 (5) = 3.10-3^. 

One can notice on Fig. [5] (left) that the graphs of these two functions intersect away from zero. This im- 
plies that depending on the dilution rate D 7 the corresponding value of Ai can be less or greater than X 2 . 
The input concentration Si n has been chosen equal to 20 and the initial state vector has been kept equal to 
(^(O),^!^),^^)) = (5,2,3). A simple calculation shows that for Q = 2.1CT 4 , one has A 1 = A 2 = f < S m . 
Then for this choice of fii, [i 2 and Q the dynamical system ^ admits positive non- hyperbolic equilibriums. 
We aim to study the numerical evolution of the species concentrations about those equilibriums in both 
models, depending on the choice of the flow Q. For this purpose, we have plotted the graphs of the species 
concentrations at steady state given by MIN3P and the chemostat model in the input control volume (see 
Fig. [HJ right), denoting by B* Mi (resp. BJ.i) f° r species i = 1,2. We observe that for Q < 10 -5 , both 
simulations present almost the same solutions. When Q increases a difference between the models begins 
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to appear until we detect a wrong prediction of the species that avoids the extinction. This happens for 
15.HT 5 < Q < 4.1CT 4 . For Q > 4.10" 4 , both simulations show again almost the same solutions. 

So we have observed another significant difference between simulations of reactive transport and chemo- 
stat models about a bifurcation point, where the species that avoid extinction switches. 




Figure 5: Two species in competition in a single chemostat 



To study the effect of a spatial discretization in the case of two species, we have chosen a specific example 
of twenty perfectly-mixed tanks that are interconnected in series. The volumes Vi(i — 1, ...,20), the input 
concentration Si n , the flux Q and the specific growth rates /ii(.) and ^(O are chosen as follows, in such 
a way that the species B\ passes from the wash-out state (in V\, V% and V3) to a coexistence state (in 
^,1 = 4, ...,20). 

Vi = to.io- 2 s 

y 2 = 9.10- , l( ^) = 4.629.10- — 

V 3 = 1.10- 2 m{s) = 6.944.10-5^^ 



v 4 = 11.10- 2 

V 5 = 9T0- 2 

V 6 = --- = V 20 = 4.10- 



18 + S 
Q = 0.3587.10- 5 

S in = 19.25 



Under MIN3P we have discretized the domain into twenty control volumes, and have compared the 
solutions of both models in each control volume at steady state. 

On Figure [6] (right), one can observe that the substrate concentrations computed at steady state with 
both models are almost identical in the first control volumes, up to the third one where the solutions are 
different, and become radically divergent in further control volumes. Here, the solution computed by MIN3P 
does not predict the coexistence of two species... Let us underline that the species coexistence is no longer 
a pathological case when chemostats of different volumes are connected in series (see for instance [4U1 135j ). 

To summarize, we have shown that under certain circumstances the issue of the winner of the competition 
between the two species are predicted radically different by the simulations of both models. 
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S N th control volume 



Figure 6: Two species in competition in a series of chemostats 



4 Conclusion 

Our results show a possible inaccuracy of numerical reactive transport models in the simulation of the 
dynamics of simple ecosystem of chemostat type. Our objective was not to challenge the model MIN3P for 
its ability to simulate complex problems involving mass flow and multicomponents reaction networks, but 
we raise the fact that the numerical accuracy of the model MIN3P depended on the evolution of microbial 
species. When the hydrodynamic conditions make the system close to a washing-out of one or several 
microbial species, or to the coexistence of species, we observe numerical bias in the computation of the 
solution that leads to radically different predictions. Consequently one may wonder if the numerical issue 
found here in simple systems (a single substrate and one or two microbial species, advective transport) prevails 
in more complex systems, when a network of kinetically-controlled reactions is considered to simulate for 
instance remediation problems in ground- waters. We believe that a study of the eigenvalues of the linearized 
dynamics about steady states is important, for detecting possible numerical inaccuracy. 



Appendix 

For the study of the behavior of autonomous dynamical systems in R", 

f-W, (3) 

with F £ C (R n ), one usually determine first its equilibrium points (denoted X*) as solutions of f(X) = 0, 
and then study the eigenvalues of the linear dynamics 

called the linearization of ^ about X*, where J(X*) is the Jacobian matrix of F at X*. If all the eigenvalues 
Vi (i — 1, ...n ) of J(X*) have nonzero real parts, then we said that X* is hyperbolic. When at least one of 
its eigenvalues has a zero real part, then we said that X* is non- hyperbolic (see for instance the reference |16|). 
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We recall now the usual definitions of stability and a main result allowing to conclude about the nature 
of an equilibrium (see the reference |16| for more details). 



Definition 4.1 The equilibrium X* is said to be stable if for every e > there exists S > such that 

\\X(t ) - X*\\ < <5 ^ \\X(t) - X*\\ < e Vi > t Q . 
If this condition its not satisfied, the equilibrium is unstable. 

Definition 4.2 The equilibrium X* is said to be (locally) exponentially stable if for every e > there exists 
three real numbers a > 0, b > and 5 > such that 

\\X{t ) - X*\\ <6^ \\X(t) - X*\\ < a \\X(t Q ) - X*\\ e~ bt Vi > t . 

Theorem 4.3 

• If all the eigenvalues of J(X*) have negative real parts then the equilibrium point X* is exponentially 
stable. 

• If one of the eigenvalues of J(X*) has a positive real part then the equilibrium point X* is unstable. 

Remark. If the Jacobian matrix J(X*) has at least one eigenvalue with zero real part, then we need to use 
other results to conclude about the behavior of the trajectory of the system pi). 

For the chemostat model 0, one has the following result. 
Proposition 4.4 Denote by A the solution when it exists of fx(S) = D. 

• If D < fi(Si rl ), the system £7J) admits two equilibrium points given by Ei(Si n , 0), which is unstable and 
E2(X,k(Si n — X)), which is locally exponentially stable. 

• If D > fj,(Si n ), the only non-negative equilibrium point it is E\ which is locally exponentially stable, 
excepted for the case D = fi(Si n ) for which it is non-hyperbolic. 

Proof. We can easily verify that for any t > 0, the trajectories of remains in the first positive orthant, 
and are bounded: one can straightforwardly write 

^ + k^ = D{kS m ~B-kS) 
dt dt 

from which one deduces that 1 1— > B{t) + kS(t) is bounded that the trajectories of the system are bounded. 
Determining the equilibrium points of the system ([I]) amounts to solve the following system 

.^ B + D(S m -S) = 

k (4) 

(/i{S) - D)B = 0. 

The wash-out equilibrium point Ei(Si n , 0) is always solution, and there is a possibility of another equilibrium 
point E2(S*, k(Si n — S*)) with S* = A, when S* < Si n . To study the stability of these equilibrium points, 
we write the Jacobian matrix of the system ([lj : 



J(S,B) = 



( AS) B D MCg) 
k k 

\ fi>(S)B fi(S) - D 
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whose eigenvalues are 



and 



vi(S,B) — -D < 
v 2 (S,B) = »(S)-D-^p-B. 



At Ei, one has v 2 (Ei) = fi(Si n ) — D and at the non-trivial equilibrium _E 2 (when it exists), one has v 2 (E 2 ) = 

K 

So, when fJ,(Si n ) < D, we conclude that the non-trivial equilibrium does not exist and we obtain that E\ is 
locally exponentially stable. When fi(Si n ) > D, E% is unstable and E 2 is locally exponentially stable. For the 
particular case (J,(Si n ) = D, the non-trivial equilibrium does not exist and E\ is a non-hyperbolic equilibrium. 

For the chemostat model ^ with two species, one has the following result. 

Proposition 4.5 Denote by Aj the solution (when it exists) of fJ,i(Xi) — D. Under the condition that Si n > 
max(Ai, Aa), the system |Ip admits three equilibrium points given by Eo(Si n ,0, 0), £i(Ai, k\{Si n — Ai), 0) and 
^2(^2, 0, k2(Si„ — A2)). Furthermore, one has 

• when Xi < Xj, Ei is locally exponentially stable and Eq and Ej are both unstable. 

• when Ai = A2 then E\ = E 2 is a non-hyperbolic equilibrium point. 
Proof. As before, the equilibrium points are given by the following system 

^ s ) Bl _f^l B2 + D{Sin -s) =0 



fci k 2 
{fi 1 {S)-D)B 1 

(MS)~D)B 2 



= 
= 



(5) 



One can find that there exist at most three equilibrium points Eo(Si n , 0, 0), E\(Xi, k\{Si n — Ai),0) and 
E 2 (X 2 ,0,k 2 (S in -X 2 )). 

For the study of their stability, we write for convenience the dynamics in variables [Z, B\,B 2 ) with 



dZ 
~dt 
dB t 
~dT 
dB 2 
dt 



Bi 
Ki ' 


B 2 
V K 2 


+ s . 


D(S in 


-z) 






-Bi 
Ki 


B 2 
K 2 




Si 


B 2 
K 2 



) - D)B 1 
) - D)B 2 



In these coordinates, the Jacobian matrix takes the following form: 

l-B 

* -^l Bl +^{S)-D 



J{Z, B U B 2 ) 



k 2 



Bi 



V 



/4(g) 
fci J 



.^B 1 + f , 2 (S)-D ) 



At Eq, we can check that J{Eq) admits three eigenvalues: 

v 1 (E Q ) = -D <0, v 2 (E ) = m(S in )-D > 0, v 3 (E ) = n 2 (S in ) - D > 
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Thus E is unstable. At E\, the eigenvalues are 



Vl (£1 ) = -£>< 0, 



v 2 (E 1 ) 



A*i(Ai) 



fci(Sj. 



Ai)<0, v 3 (£i) = m 2 (Ai) - D 



and symmetrically for E 2 : 



vi(E 2 ) = -D < 0, 



V 2 (E 2 ) 



/4(A 2 ) 

fc 2 



MS, 



A 2 )<0, v 3 (£ 2 ) =/ii(A 2 )-D ■ 



One can notice that Aj < Aj <^ v 3 (E i - ) < and conclude that when Aj < Aj, is locally exponentially 
stable and Ej is unstable. 

For the particular case of Ai = A 2 we find that v 3 (Ei) = v 3 (E 2 ) — 0. Therefore in this case E\ (that 
coincides with E 2 ) is a non- hyperbolic equilibrium point. 
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